home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_3.4 / peak / convsqrs.c next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  7.2 KB  |  206 lines
  1. /* 
  2.  * convsqrs.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /* CONVSQRS:    function performs converging squares algorithm,
  10.  *                            which finds a peak by iteratively choosing highest
  11.  *                              density peak from low to higher resolutions
  12.  */
  13.  
  14. #define LT_GT >                 /* greater than finds upward peaks for PEAKUP */
  15. #define ELEMENT(P,N,I,J) (( (int) (*(P + I * N + J))) & 0x0ff)
  16.     /* returns int from byte array at address P, coord.s I,J
  17.      * from an array of length N */
  18. /* macro finds sum of elements in row of image */
  19. #define SUMROW(IMAGE,WIDTH,ROW,SCOL,ECOL,SUM) \
  20.     SUM = 0; \
  21.     pImage = IMAGE + (ROW) * WIDTH + SCOL; \
  22.     for (col = SCOL; col <= ECOL; col++, pImage++) { \
  23.         SUM += ((long) *pImage) & 0x0ff; \
  24.     }
  25.  
  26. /* macro finds sum of elements in column of image */
  27. #define SUMCOL(IMAGE,WIDTH,COL,SROW,EROW,SUM) \
  28.     SUM = 0; \
  29.     pImage = IMAGE + (SROW) * WIDTH + COL; \
  30.     for (row = SROW; row <= EROW; row++, pImage += WIDTH){ \
  31.         SUM += ((long) *pImage) &0x0ff; \
  32.     }
  33.  
  34. /* macro finds subsquare of peak sum. the order of testing of subsquares
  35.  * is revolved so as to find centroid of equi-intensity plateau */
  36. #define PEAKSQR(SQRSUM,POSITION,PEAK,I,TEMP) \
  37.     TEMP = I[0]; \
  38.     I[0] = I[1]; \
  39.     I[1] = I[2]; \
  40.     I[2] = I[3]; \
  41.     I[3] = TEMP; \
  42.     PEAK = SQRSUM[I[0]]; \
  43.     POSITION = I[0]; \
  44.     if  (SQRSUM[I[1]] LT_GT PEAK){ \
  45.         PEAK = SQRSUM[I[1]]; \
  46.         POSITION = I[1]; \
  47.     }\
  48.     if  (SQRSUM[I[2]] LT_GT PEAK){ \
  49.         PEAK = SQRSUM[I[2]]; \
  50.         POSITION = I[2]; \
  51.     } \
  52.     if (SQRSUM[I[3]] LT_GT PEAK){ \
  53.         PEAK = SQRSUM[I[3]]; \
  54.         POSITION = I[3]; \
  55.     }
  56.  
  57. long
  58. convsqrs (image, width, height, rowInitial, colInitial, nInitial,
  59.           nFinal, rowPeak, colPeak)
  60.      unsigned char *image;      /* pointer to image array */
  61.      long width, height;        /* image sidelengths */
  62.      int rowInitial;            /* initial start row and column of conv. sqrs wrt image */
  63.      int colInitial;
  64.      int nInitial;              /* initial size of converging squares */
  65.      int nFinal;                /* sidelength of final converging square */
  66.      int *rowPeak;              /* row coord. of peak within final square */
  67.      int *colPeak;              /* col coord. of peak within final square */
  68. {
  69.   unsigned char *pImage;        /* pointer to image array */
  70.   long sqrSum[4];               /* sum of elements of squares 1-4 */
  71.   long col1Sum;                 /* sum of left k-2 length col.s of square */
  72.   long col2Sum;                 /* sum of right ... */
  73.   long row1Sum, row2Sum;        /* sum of top, bottom ... */
  74.   long peak;                    /* peak value in final square */
  75.   long value;                   /* long value from byte array */
  76.   int position;                 /* no. (0-3) of max. square */
  77.   int nSquare;                  /* incrementor = current sidelength reduced square */
  78.   int sqrIndex[4];              /* revolving indices of 4 subsquares */
  79.   int temp;
  80.   register int sRow, eRow;      /* start and end rows cols of output reduced sqr */
  81.   register int sCol, eCol;
  82.   register int row, col;        /* row and column incrementors */
  83.  
  84. /* initialize */
  85.   if (nFinal < 2)
  86.     nFinal = 2;                 /* square size of reduced array */
  87.   sRow = rowInitial;
  88.   sCol = colInitial;
  89.   if ((sRow + nInitial) >= height) {
  90.     nInitial = height - 1 - sRow;
  91.   }
  92.   if ((sCol + nInitial) >= width) {
  93.     nInitial = width - 1 - sCol;
  94.   }
  95.   eRow = sRow + nInitial;
  96.   eCol = sCol + nInitial;
  97.  
  98. /* initialize revolving square indices */
  99.   sqrIndex[0] = 0;
  100.   sqrIndex[1] = 1;
  101.   sqrIndex[2] = 2;
  102.   sqrIndex[3] = 3;
  103.  
  104. /* reduce initial (nImage x nImage square to (nImage-1)x(nImage-1) square */
  105.  
  106. /* find sums of start and end columns and rows */
  107.   SUMCOL (image, width, sCol, sRow + 1, eRow - 1, col1Sum);
  108.   SUMCOL (image, width, eCol, sRow + 1, eRow - 1, col2Sum);
  109.   SUMROW (image, width, sRow, sCol + 1, eCol - 1, row1Sum);
  110.   SUMROW (image, width, eRow, sCol + 1, eCol - 1, row2Sum);
  111.  
  112. /* find sums of the 4 (nImage -1)x(nImage-1) squares, and peak square */
  113.   sqrSum[0] = col1Sum + row1Sum + ELEMENT (image, width, sRow, sCol);
  114.   sqrSum[1] = row1Sum + col2Sum + ELEMENT (image, width, sRow, eCol);
  115.   sqrSum[2] = col1Sum + row2Sum + ELEMENT (image, width, eRow, sCol);
  116.   sqrSum[3] = row2Sum + col2Sum + ELEMENT (image, width, eRow, eCol);
  117.  
  118.   PEAKSQR (sqrSum, position, peak, sqrIndex, temp);
  119.  
  120. /* reduce from (nImage-1)x(nImage-1) square to (nFinal x nFinal) square
  121.  * sequentially, for squares of peak density */
  122.  
  123.   for (nSquare = nInitial - 1; nSquare > nFinal; nSquare--) {
  124.     /* calculate end row and column sums for peak square */
  125.     switch (position) {
  126.     case 0:                    /* top left square */
  127.       eRow -= 1;
  128.       eCol -= 1;
  129.       col1Sum = col1Sum - ELEMENT (image, width, eRow, sCol);
  130.       SUMCOL (image, width, eCol, sRow + 1, eRow - 1, col2Sum);
  131.       row1Sum = row1Sum - ELEMENT (image, width, sRow, eCol);
  132.       SUMROW (image, width, eRow, sCol + 1, eCol - 1, row2Sum);
  133.       break;
  134.     case 1:                    /* top right square */
  135.       sCol += 1;
  136.       eRow -= 1;
  137.       SUMCOL (image, width, sCol, sRow + 1, eRow - 1, col1Sum);
  138.       col2Sum = col2Sum - ELEMENT (image, width, eRow, eCol);
  139.       row1Sum = row1Sum - ELEMENT (image, width, sRow, sCol);
  140.       SUMROW (image, width, eRow, sCol + 1, eCol - 1, row2Sum);
  141.       break;
  142.     case 2:                    /* bottom left square */
  143.       sRow += 1;
  144.       eCol -= 1;
  145.       col1Sum = col1Sum - ELEMENT (image, width, sRow, sCol);
  146.       SUMCOL (image, width, eCol, sRow + 1, eRow - 1, col2Sum);
  147.       SUMROW (image, width, sRow, sCol + 1, eCol - 1, row1Sum);
  148.       row2Sum = row2Sum - ELEMENT (image, width, eRow, eCol);
  149.       break;
  150.     case 3:                    /* bottom right square */
  151.       sRow += 1;
  152.       sCol += 1;
  153.       SUMCOL (image, width, sCol, sRow + 1, eRow - 1, col1Sum);
  154.       col2Sum = col2Sum - ELEMENT (image, width, sRow, eCol);
  155.       SUMROW (image, width, sRow, sCol + 1, eCol - 1, row1Sum);
  156.       row2Sum = row2Sum - ELEMENT (image, width, eRow, sCol);
  157.       break;
  158.     }
  159.  
  160. /* find sums of 4 (nSquare-1)x(nSquare-1) squares and peak square */
  161.     sqrSum[0] = col1Sum + row1Sum + ELEMENT (image, width, sRow, sCol);
  162.     sqrSum[1] = row1Sum + col2Sum + ELEMENT (image, width, sRow, eCol);
  163.     sqrSum[2] = col2Sum + row2Sum + ELEMENT (image, width, eRow, sCol);
  164.     sqrSum[3] = row2Sum + col2Sum + ELEMENT (image, width, eRow, eCol);
  165.  
  166.     PEAKSQR (sqrSum, position, peak, sqrIndex, temp);
  167.   }
  168.  
  169. /* find coordinates of (nFinal x nFinal) square */
  170.   switch (position) {
  171.   case 0:                      /* top left square */
  172.     eRow -= 1;
  173.     eCol -= 1;
  174.     break;
  175.   case 1:                      /* top right square */
  176.     sCol += 1;
  177.     eRow -= 1;
  178.     break;
  179.   case 2:                      /* bottom left square */
  180.     sRow += 1;
  181.     eCol -= 1;
  182.     break;
  183.   case 3:                      /* bottom right square */
  184.     sRow += 1;
  185.     sCol += 1;
  186.     break;
  187.   }
  188.  
  189. /* find peak element and its position within final square */
  190.   peak = ELEMENT (image, width, sRow, sCol);
  191.   *rowPeak = sRow;
  192.   *colPeak = sCol;
  193.  
  194.   for (row = sRow; row <= eRow; row++) {
  195.     for (col = sCol; col <= eCol; col++) {
  196.       if ((value = ELEMENT (image, width, row, col)) LT_GT peak) {
  197.         peak = value;
  198.         *rowPeak = row;
  199.         *colPeak = col;
  200.       }
  201.     }
  202.   }
  203.  
  204.   return (peak);
  205. }
  206.